home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / answrbok / 2_2.lha / 2_2 / rsqrt.c < prev   
Text File  |  1993-08-08  |  2KB  |  69 lines
  1. * Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
  2. * The C++ Answer Book */
  3. * Tony Hansen */
  4. * All rights reserved. */
  5. *ident    "@(#)cfront:lib/complex/sqrt.c    1.4" */
  6.  include "complex.h"
  7.  
  8.  
  9. define    SQRT_DANGER    1e17
  10.  define PERIL(t)    (t > SQRT_DANGER || (t < 1/SQRT_DANGER && t != 0))
  11.  
  12. *
  13. *
  14. * 7-25-83, note from Leonie Rose -
  15. * Stu Feldman says that the peril tests for the following function
  16. * are "acceptable" for now, but certain things like
  17. * sqrt(1e10 + 1e-30*i) will cause floating exceptions.
  18. *
  19. */
  20.  
  21. omplex
  22. qrt(complex z)
  23.  
  24. complex  answer(1,4.5);
  25. double  r_old,  partial;
  26.  
  27. *
  28. Check for possible overflow, and fixup if necessary.
  29. /
  30.  
  31. double x = abs(z.re);
  32. double y = abs(z.im);
  33.  
  34. if (x > y && PERIL(x)) {
  35.             z.im /= x;
  36.            z.re /= x;  /* z.re is replaced by 1 or -1 */
  37.             partial = sqrt(x);
  38.         }
  39. else if PERIL(y) { 
  40.             z.im /= y;  /* roles of z.re, z.im reversed from previous */
  41.             z.re /= y;
  42.             partial = sqrt(y);
  43.         }
  44. else partial = 1;
  45.  
  46. *
  47. Main computation:
  48. Use half angle formulas to compute angular part of the square root.
  49. The sign of sin_old is the same as that for sin_new, which means that the
  50. upper half plane gets mapped to the first quadrant, and
  51. the lower half plane to the fourth quandrant.
  52. /
  53.  
  54.  
  55. if (r_old = sqrt(z.re*z.re + z.im*z.im)) {
  56.            double r_new = partial * sqrt(r_old);
  57.  
  58.            double cos_old = z.re/r_old;
  59.            double sin_old = z.im/r_old;
  60.            double cos_new = sqrt( (1 + cos_old)/2 );
  61.           double sin_new = (cos_new == 0)? 1 : sin_old/(2*cos_new);
  62.  
  63.            answer.re = r_new * cos_new;
  64.            answer.im = r_new * sin_new;
  65.        }
  66.  
  67. return answer;
  68.  
  69.